disp('fmdlf');

an = @(x, n) 1/2/pi*n./((n/2).^2 + x.^2);

vx = linspace(0,10,400);
m = 4;
vf = an(vx, m);
n = 2;
vx_out = vx;

[vf_out, vf_out0] = fn_anf(n, vx, vf, vx_out);

close(figure(1));
figure(1);
plot(fn_ex1(vx), fn_ex2(vf_out), '-ro');
hold on;
plot(fn_ex1(vx_out), fn_ex2(an(vx_out, m+n)), '-b*');

diss = 1/sqrt(length(vx_out))*norm(vf_out - an(vx_out, m+n));
fprintf('error = %e \n',  diss);

diss = 1/sqrt(length(vx_out))*norm(vf_out0 - an(vx_out, m+n));
fprintf('error0 = %e \n',  diss);

vep = vx.^2  - 1;
F = @(x,T) -T*log(1+exp(-x/T));
vf = F(vep, 0.2);
close(figure(2));
figure(2);
plot(fn_ex1(vx), fn_ex2(vf), '-ro');
[vf_out, vf_out0] = fn_anf(n, vx, vf, vx_out);
figure(2); hold on;
plot(fn_ex1(vx_out), fn_ex2(vf_out), '-bo');

diss = 1/sqrt(length(vx_out))*norm(vf_out - an(vx_out, m+n));
fprintf('error = %e \n',  diss);

diss = 1/sqrt(length(vx_out))*norm(vf_out0 - an(vx_out, m+n));
fprintf('error0 = %e \n',  diss);


vep = vx.^2  - 1;
F = @(x,T) -T*log(1+exp(-x/T));
T = 0.2;
vf = F(vep, T);
vu = -T*fn_anf(2, vx, vf, vx_out);
close(figure(4));
figure(4);
plot(fn_ex1(vx_out), fn_ex2(vu), '-bo');
a = fn_inveta(2, 1);
tv = log((1 + exp(-vu)*a)./(1 + a));
figure(4);
hold on;
plot(fn_ex1(vx_out), fn_ex2(tv), '-r*');

n = 2;
m = 4;
vu_out = fn_Tnm(n, m, fn_inveta(m,1), vx, vu, vx_out);
figure(4);
plot(fn_ex1(vx_out), fn_ex2(vu_out), '-b*');
disp(fn_inveta(m,1));

figure(44);
plot(fn_ex1(vx_out), fn_ex2(vu_out), '-b*');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function vu_out = fn_Tnm(n, m, a, vx, vu, vx_out)

vu_out = 0; 

vf = log((1 + exp(-vu)*a)./(1 + a));

if m == n
    
    for i = 2:2:(2*n-2)
        
        vf_out = fn_anf(i, vx, vf, vx_out);
        vu_out = vu_out + 2*vf_out;
        
    end
    
    vu_out = vu_out + fn_anf(2*n, vx, vf, vx_out);
    
else
    
    vu_out = vu_out + fn_anf(abs(m-n), vx, vf, vx_out);
    vu_out = vu_out + fn_anf(m+n, vx, vf, vx_out);
    
    for i = (abs(m-n)+2):2:(m+n-2)
        
       vu_out = vu_out + 2*fn_anf(i, vx, vf, vx_out); 
        
    end
    
end

end




function res = fn_inveta(m, h0)

res = sinh(h0/2)^2/(sinh(h0/2*(m+1))^2 - sinh(h0/2)^2);

end



function [vf_out, vf_out0] = fn_anf(n, vx, vf, vx_out)

an = @(x) 1/2/pi*n./((n/2).^2 + x.^2);

vk = vx_out;
vp = vx;

mk = transpose(vk)*ones(1,length(vp));
mp = ones(length(vk),1)*vp;
ma = an(mk - mp) + an(mk + mp);

vpd = vp(2:end) - vp(1:end-1);
vp_nt = zeros(1,length(vp));
vp_nt(1:end-1) = 1/2*vpd;
vp_nt(2:end) = vp_nt(2:end) + 1/2*vpd;

vf_out = ma*diag(vp_nt)*transpose(vf);
vf_out = transpose(vf_out);
vf_out0 = vf_out;

xc = vx(end);
alpha = xc^2*vf(end);
vf_out_addon = alpha ...
    *(fn_int_tail(vk, xc, n) + fn_int_tail(-vk, xc, n));

vf_out = vf_out + vf_out_addon;

end


function vx_ex1 = fn_ex1(vx)

vx_ex1 = [-vx(end:-1:2), vx];

end

function vx_ex2 = fn_ex2(vx)

vx_ex2 = [vx(end:-1:2), vx];

end


function res = fn_int_tail(k, kc, n)

res = 8/pi*n*log(4)*k./(4*k.^2 + n^2).^2 ...
    + (4*k.^2 - n^2)*2./(4*k.^2 + n^2).^2 ...
    + 2/pi*n./(4*k.^2 + n^2)./kc ...
    - 4/pi*(4*k.^2 - n^2)./(4*k.^2 + n^2).^2.*atan(2*(k + kc)./n) ...
    + 8/pi*n*k./(4*k.^2 + n^2).^2.*log(kc.^2./(n^2 + 4*(k + kc).^2));

end
